function s_functions(prob, lambda, D0, chi, phi, R1, delta, dist_s, s_min, s_max, dist_d, d_min,d_max, R1_grid, delta_grid, optfig)
% Illustration of s_hat and s_star
R1_saved = R1;

% Varying R1

s_hatR  = zeros(size(R1_grid));
s_starR = zeros(size(R1_grid));
for i=1:length(R1_grid)
    R1_i       = R1_grid(i);
    
    s_hatR(i)  = fn_s_hat(R1_i,lambda,s_min,s_max,phi);
    s_starR(i) = fn_s_star(delta,R1_i,D0,lambda,s_min,s_max,dist_d,d_min,d_max,phi);
end

prob_fail_R   = dist_s.cdf(s_hatR) + prob*(dist_s.cdf(s_starR) - dist_s.cdf(s_hatR)); % Failure probability
prob_fail_F_R = dist_s.cdf(s_hatR); % Fundamental failure

%% Varying delta

R1 = R1_saved;
s_hat    = zeros(size(delta_grid));
s_star   = zeros(size(delta_grid));
for i=1:length(delta_grid)
    del_i     = delta_grid(i);
    
    s_hat(i)  = fn_s_hat(R1,lambda,s_min,s_max,phi);
    s_star(i) = fn_s_star(del_i,R1,D0,lambda,s_min,s_max,dist_d,d_min,d_max,phi);
end

prob_fail   = dist_s.cdf(s_hat) + prob*(dist_s.cdf(s_star) - dist_s.cdf(s_hat)); % Failure probability
prob_fail_F = dist_s.cdf(s_hat); % Fundamental failure

fig_s_functions(R1_grid, s_hatR, s_starR, delta_grid, s_hat, s_star, s_min, s_max, ...
    prob_fail, prob_fail_F, prob_fail_R, prob_fail_F_R , optfig)

%% Varying s

s_grid = (s_min:0.001:s_max)';

tax    = zeros(size(s_grid));
alphaF = zeros(size(s_grid));
alphaN = zeros(size(s_grid));
for i=1:length(s_grid)
    s = s_grid(i);
    
    tax(i)    = fn_tax(delta,R1,D0,chi,s,phi,dist_d,d_min,d_max);
    alphaF(i) = fn_alphaF(delta,R1,D0,chi,s,phi,dist_d,d_min,d_max);
    alphaN(i) = fn_alphaN(R1,lambda,phi,s);
end

fig_s_paper(s_grid, tax, alphaN, alphaF, optfig)

end
